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' (DPM) and decentralized Lanczos algorithm (DLA) - for distributed computation of one (the largest) or 

I multiple eigenvalues of a sample covariance matrix over a wireless network. The proposed algorithms, 

based on sequential average consensus steps for computations of matrix- vector products and inner vector 
products, are first shown to be equivalent to their centralized counterparts in the case of exact distributed 
CO ' consensus. Then, closed-form expressions of the error introduced by non-ideal consensus are derived 
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Abstract 

In this paper we derive and analyze two algorithms - referred to as decentraUzed power method 



for both algorithms. The error of the DPM is shown to vanish asymptotically under given conditions 
on the sequence of consensus errors. Finally, we consider applications to spectrum sensing in cognitive 

m 

i radio networks, and we show that virtually all eigenvalue-based tests proposed in the literature can be 



implemented in a distributed setting using either the DPM or the DLA. Simulation results are presented 
that validate the effectiveness of the proposed algorithms in conditions of practical interest (large-scale 



I networks, small number of samples, and limited number of iterations). 
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I. Motivation and Related Work 

Computing the eigenvalues of sample covariance matrices is a fundamental problem in signal 
processing, with applications including multi- sensor spectrum sensing in cognitive radio networks 
(see Section VI). Given the increasing popularity of dense, large-scale wireless sensor networks, 
applications of eigenvalue-based inference techniques in distributed settings are of great interest. 
However, most of the eigenvalue-based techniques proposed in the existing literature assume 
a centralized architecture, where the samples received by different nodes are forwarded to a 
fusion center that is in charge of constructing the sample covariance matrix and computing the 
relevant test statistics. This traditional architecture has several drawbacks: it requires a fusion 
center with high computational capabilities, therefore it does not support applications in which 
different nodes may be chosen as fusion centers at different times; relying on one node only, 
it is vulnerable to hardware failures, Byzantine faults, or attacks from malicious users; it is not 
efficient in case of multi-hop networks, where some nodes may be many hops away from the 
fusion center; and it lacks scalability, because a growing number of nodes in the network may 
result in congestion of the communication channel with the fusion center. 

For these reasons, we seek a decentralized method to compute the eigenvalues of sample 
covariance matrices over a wireless network, such that the computational effort is distributed 
across multiple nodes and the many-to-one communication protocol is replaced by a more scal- 
able neighbor-to-neighbor protocol. In this paper we propose two solutions based on decentralized 
implementations of iterative eigenvalue algorithms - the power method (PM, see Section III-B) 
and the Lanczos algorithm (LA, see Section III-C). Decentralized versions of such algorithms 
are obtained by applying average consensus (AC, see Section III-A) as a subroutine to perform 
those computations that involve combining the data of different nodes. Once local estimates of 
the eigenvalues of interest are computed at every node, statistical tests for signal detection can 
be performed locally by each node. 

The contribution of this paper is related to that of [3], where a decentralized algorithm based 
on the Oja-Karhunen recursion is proposed to track the eigenvectors of a covariance matrix. 
Our work adopts the same idea of computing inner vector products through AC (and further 
extends this approach to matrix-vector products), but differs from [3] in the following sense: 
(i) we compute eigenvalues instead of (possibly, in addition to) eigenvectors, thus adapting the 
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well-established theory of eigenvalue-based detection to decentralized network settings; (ii) we 
focus on detection (decision based on N received samples per sensor) instead of sequential 
tracking (update of eigenvector estimates at every new sample). This makes our approach 
suitable for spectrum sensing and other signal detection applications. A similar methodology for 
distributed matrix multiplication via AC has been recently used also in [4], where a decentralized 
expectation-maximization algorithm is derived for static linear Gaussian models. 

Other recent related works are [5] and [6], both proposing decentralized implementations of 
principal component analysis (PCA) over wireless networks. The first approach [5] relies on 
the assumption of decomposable Gaussian graphical models, i.e., it requires data sample to 
be (i) multivariate Gaussian distributed, and (ii) decomposable into two or more conditionally 
independent "cliques". The global eigenvalue decomposition (EVD) problem is thus broken down 
into a sequence of local (clique- wise) EVD subproblems. The second approach [6] combines the 
power method with the concept of sparsification to achieve an efficient distributed computation 
of the eigenvectors and eigenvalues of a symmetric matrix. However, this approach is based 
on the assumptions that (i) each node has access to a full row of the matrix (which is not the 
case with the sample covariance matrix considered in our model, see Section II), and (ii) the 
network graph is completely connected (i.e., a direct link exists between any pair of nodes). Our 
approach, on the contrary, does not require any of the aforementioned assumptions. 

The paper is organized as follows. A formal statement of the problem is provided in Section 
II. Section III contains mathematical preliminaries about the algorithms used in this work (AC, 
PM, and LA). We then introduce the proposed algorithms in Section IV and analyze their 
performance and complexity in Section V. We finally discuss two practical applications of the 
proposed algorithms in Section VI and present numerical results in Section VII. Concluding 
remarks are provided in Section VIII. 

II. Problem Statement 

Consider a wireless network consisting of K sensor nodes. During a given time interval 
(sensing period), each node collects complex signal samples. The global received sample 
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matrix is denoted by 



Y=[y(l),...,y(iV)] 



^KxN 



(1) 



where symbols y(-) G C and y[-] G C are used to denote, respectively, the columns and 
(transpose) rows of Y. Physically, column y{n) contains the samples received by all nodes at 
time n, while row y[A;]^ contains all samples available at node k at the end of the sensing period. 
We then define the sample co variance matrix as 



R^-YY 



H 



(2) 



Let Ai > ... > \k > be the eigenvalues of R, without loss of generality sorted in 
decreasing order, and ui,...,ux the corresponding eigenvectors. The problem addressed in 
this work can be stated as follows: how can a network compute (or estimate) one or more of 
the above eigenvalues without a fusion center that collects all samples Y, and without explicitly 
constructing the sample covariance matrix R? Before presenting the proposed solution, we 
introduce some preliminary definitions and basic concepts about distributed consensus, power 
method, and Lanczos algorithm. The notation used throughout the paper is summarized in Table 
I. 



in. Preliminaries 

A. Average Consensus 

Assume that the network nodes and their links form a connected undirected graph. Under 
such an assumption, it is possible to define a distributed AC algorithm over the network. By 
distributed AC we mean any algorithm whose output for all nodes converges to the average of the 
initial values of the individual nodes. A large variety of AC algorithms have been proposed in the 
literature (see [7] for a survey), both with synchronous [8], [9] and asynchronous protocols [10]- 
[12]. Extensions to noisy message exchange and link failures include [13]-[17], and methods 
for consensus acceleration have been proposed in [18]-[20]. It is worth noting that, under the 
assumption of fixed network topology and noiseless communication, exact consensus can be 
achieved in a finite number of steps: see for example [21]. 
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Symbol 


Definition 






l|v|| 


Euclidean norm of vector v 






T * H 


Transpose, complex conjugate, and Hermitian operators 






a{n), a[k]'^ 


Respectively, column n and row k of matrix A 






0, 


Element- wise vector multiplication and division 








Identity matrix and column vector of ones of size m (subscript is omitted when not ambiguous) 


diag(A) 


Column vector with diagonal elements of a square matrix A 


: i.e., diag(A) = (A I) • 


1 


Diag(v) 


Square matrix with vector v as main diagonal: i.e., Diag(v) 


= (v-l^)0l 




RN(A) 


Vector of l'^ square norms of the rows of A: i.e., RN(A) = 


[||a[l]f,...,||a[Alf]^ = 


(A* A) • 1 


aC(-) 


AC function with global input/output for all nodes (Eq. 4) 








m = input vector size, t = AC iterations or running time (no superscript = ideal AC) 




AC^[fc](-) 


AC function with local input/output for node k (Eq. 6) 







TABLE I 
Notation. 



In this paper we do not adopt a specific AC algorithm, but rather take a general approach. We 
model the result of a generic AC routine by a function AC^^ : C^^"* — > C^^™, where m is the 
size of the input vectors at each node^ and t is the number of iterations or the averaging time.^ 
Denote the global input (or initial value) matrix by 

Zo = [zo(l),...,Zo(m)] = ; 

_ MKf 

defined in analogy with (1), such that the k-th row represents the samples available at node k 
Then, the output of the AC function at time t is defined by 



(3) 



AC^,(Zo)^-ll'Zo + E 



(4) 



'An AC algorithm with vector inputs involves exchanging m scalar numbers at every iteration, and returns the average of the 
input vectors over all network nodes. As such, it involves the same number of messages as scalar AC, but with a larger payload. 

^For sake of generality, we do not specify whether the adopted AC routine should be synchronous or asynchronous. In the tirst 
case t is a discrete number of iterations; in the second case it is the running time of the algorithm, that is related probabilistically 
to the number to clock ticks (see for example [11]). 
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where 1 is a column vector of ones of size K, and 

Ei = [et(l),...,ei(m)]GC^""^ (5) 

is an error term depending on the AC time t and on the specific AC method. In general, Ej can 
be assumed to be either bounded (at least statistically, see [9]-[12], [14]-[16], [19]) or equal to 
zero (in the case of finite-time AC [21]). 

In the following, it is sometimes convenient to express the input and the output of AC for 
a single node k. For this purpose, we define a function AC^[A;] : C^^'" — )■ C^^™ having as 
input/output the k-th rows of the input/output matrices of the global function AC^(-) defined in 
(4). Thus, we can write 

z,[A;]^ = AC*^[A;](zo[A;]^). (6) 

When AC is applied to scalar arguments (m = 1), the global input-output relation is written as 
zj = kC{{zo), with zo,zt G C^^. Finally, if Et = in (4), we call the AC routine "ideal" and 
we use the notation ACm(-) (without superscript). 



B. Power Method 

The PM is a well-known iterative algorithm that, given a square matrix, converges to the 
eigenvector associated with the largest eigenvalue of the matrix [22]. The iteration, applied to 
the sample covariance matrix R, can be written as 

vu+i) = R-Vo), (7) 

where v^o, G C^*" is an arbitrary starting vector. By (7), the j-th iteration can be written as 
V(j) = R-' V(o) and, for j oo, the vector V(j) converges to a multiple of the eigenvector ui. The 
convergence rate of the algorithm is 0((A2/Ai)^^) [22]. After M iterations, the largest eigenvalue 
of R can be approximated by 

A, = ,8) 

C. Lanczos Algorithm 

A more sophisticated eigenvalue algorithm is the LA, originally proposed by Lanczos in [23]. 
The LA is applicable only to symmetric or Hermitian matrices (which is the case of the sample 
covariance matrix R considered here), but provides estimates of multiple eigenvalues (the number 
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of estimated eigenvalues depends on the iterations of the algorithm) and has faster convergence 
than the PM [24, Chapter 6]. The advantage of the LA over the PM lies in the fact that the LA 
takes into account, at every iteration j, the complete Krylov subspace 



/Cj(R, V(o)^ 



span{V(o),RV(o), . . . ,R-'V(o)} 



(9) 



whereas the PM only considers the last term R-'Vfo). The LA has been thoroughly studied by 
Paige [25]-[27]. In particular, several computational variants are compared in [25], showing 
that some are numerically more stable than others. Among the "stabler" variants, we adopt the 
one named 'A(l,7)" in [25] or equivalently "A2" in [26]. The same version of the algorithm is 
presented in [28, p. 651]. This variant is convenient in view of the decentralized implementation 
that will be developed in Section IV-B. 

The derivation of the algorithm can be briefly outlined as follows. Due to the Hermitian 
structure of R, for a given M < K, we can write 



RV = VT 



(10) 



where the columns of V G C^^^^ are mutually orthogonal unit-norm vectors and T G C^^^^ 
is a tridiagonal matrix. If M = K, matrices T and R are similar, so their eigenvalues are the 
same. However, as Lanczos first noted, the eigenvalues of T (sometimes referred to as "Ritz 
values") turn out to be excellent approximations of the eigenvalues of R even when M < K. 
The LA is thus defined by iteratively equating the columns of RV to those of VT. If we let 



a 



(1) 



/3, 



(2) 



(2) 



a 



(2) 



/3, 



(3) 



tt(M-l) P(M) 



/3, 



(M) Cl{A/) 



(11) 



the j-th iteration of the LA can be written as 



«0) = vJ)RV(,) 



w, 



(12) 
(13) 
(14) 
(15) 
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with an arbitrary starting vector v^o) of unit norm, and = 0. The above iteration is repeated 
for j = 1 to M, thus obtaining the coefficients a^J-^ and which are necessary to construct T. 
The desired estimates Ai, . . . , \m of the eigenvalues of R are then set to be the eigenvalues of 
T. Note that the eigenvalues of a tridiagonal matrix of size M x M can be efficiently computed 
with complexity O(M^) by using the spectral bisection method [29]. 

IV. Proposed Algorithms 

We now investigate how the aforementioned PM and LA can be implemented in a distributed 
fashion over a wireless network. That is, the goal is for each node A; G {1, 2, ... , K} to compute 
local estimates {Aj[A;]} of the eigenvalue(s) of interest, namely, ? = 1 for the PM and 1 < i < M 
for the LA. Decentralized eigenvalue estimates should be as close as possible to their centralized 
counterparts: ideally, for every eigenvalue Aj of interest, we would like to have 

Xi[k] = Xi Wk. (16) 

In the following sections we show that efficient DPM and DLA schemes can be developed by 
distributing the PM and LA vector iterations - respectively, (7) and (12)-(15) - in such a way that 
the A;-th element of vector V(j), indicated by is computed by node k. This is achieved by 

iterative exchange of messages between node k and its neighbors, using AC routines. A similar 
principle was used in [3] in order to develop a distributed implementation of the Oja-Karhunen 
recursion for eigenvectors, as discussed also in Section I. Once the elements v^^-^ [k] are available 
at the K nodes, inner product and norms necessary for eigenvalue computation are performed 
again by AC, while all element-wise operations (sums, multiplications by constants) are done 
locally at each node. Next, for both PM and LA, we first rewrite the global iteration in a way 
that is amenable to decentralized computation via AC, and then we break the global iteration 
down into a sequence of algorithmic steps to be executed by individual nodes. 

A. Decentralized Power Method 

The main result for the PM vector iteration is given by the following proposition. 
Proposition 1: Given an ideal AC routine AC.(-), the PM iteration (7) can be rewritten as 

V(,^,, = ^diag{Y ■ [AC^(Diag(v*,) ■ Y)]""}. (17) 
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Proof: From (7) we can write 

Vo>i, = ^YY%,, (18) 
= ldiag(YY^V(,,l^) (19) 

= i^diag(YY^Diag(v,,,)ll'') (20) 

= i^diag[Y(ll^Diag(v*,)Y)'']. (21) 

Now, if we let Z = ill'^Diag(v* ))Y and Zq = Diag(v*))Y, it is clear from (4) that Z is the 
ideal AC output with Zq as initial matrix: 

ll^Diag(v*,,)Y = K ■ AC^(Diag(v*,,)Y). (22) 

Combining (21) with (22) yields (17). □ 
Complicated though it may look, expression (17) naturally leads to a decentralized implemen- 
tation thanks to the following properties: (i) the A;-th element of the output vector, V(j+i)\k], is 
the A;-th element of diag(YZ^), hence 

^(,+,)[A;] =z[A;]^y[A;] (23) 

which can be computed locally by node k (recall that y[k] is the vector of samples received by 
node k and z[A;] is the A^-dimensional AC output at node k)\ (ii) the input to ACiv(-) is such 
that the k-th row only contains node A;'s local data f(j)[A;]* ■ y[A;]^ (recall that has been 

computed by node k in the preceding iteration). 

Assume now that the DPM iteration (17) has been repeated M times. ^ The largest eigenvalue 
estimate Ai (8) can be computed for all nodes by two additional calls to AC, as follows from 
the following proposition. 

Proposition 2: Given an ideal AC routine AC.(-), after M iterations of (17), K local copies 
of the largest eigenvalue estimate (8) can be obtained as 

Ailx = ^RN[AC;v(Diag(v*,,)) ■ Y)] ACi(v*,,, V(„)). (24) 

where RN : C^^*" — is a function that returns the squared i'^ norms of the rows of an input 
matrix (see Tab. I). 

'Note that the number of algorithm iterations is denoted by M because M is in fact the dimension of the underlying Krylov 
subspace. This identity is more evident in the case of the LA. 
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Proof: We first write the global output for all K nodes after M iterations as 

Aili^ = (li^vJ,)RV(M)) (li^vJ,)V(A,)), (25) 
where we have used (8) and applied element-wise division. The numerator can be written as 



lv;^„Rv,,„ = llv;^,,YY^V(„, (26) 



^l|K,Y|r (27) 

ll||l^Diag(v*,jY|p (28) 

iRN[ll^Diag(v*„,)Y] (29) 

^RN[AC^(Diag(v*,jY)]. (30) 



For the denominator, we can write 



lvJ„v,M) = ll^K,,,0V(„,) (31) 
= K-ACi(v*,,,0V(,,)). (32) 

Combining (30) with (32) and simplifying K finally yields (24). □ 
Similar to (17), expression (24) can be readily implemented in a distributed manner. At the 
numerator, every node k computes an AC vector z[A;] G starting from the initial value 
^(A/)[^]* ■ y[^]^' like in the vector iteration, and takes the norm of z[A;] locally. At the 
denominator, the scalar AC input at node k is simply the local quantity |f(A/)[^]P- Element- wise 
division is then performed internally by each node. 

Thanks to the results of Propositions 1 and 2, and replacing the ideal AC routine by one with 
finite averaging time t, the DPM can be written in algorithmic form as summarized in Alg. 1. 
The averaging time or number of iterations t is assumed to be either predefined or adjusted 
online at every iteration, based on the starting vector and a target error bound. 

B. Decentralized Lanczos Algorithm 

Consider now the LA iteration (12)-(15). We note that (12) has the same structure as the 
numerator of (8), and (13) is similar to the PM iteration (7), with additional terms which can 
be computed locally provided that each node k has stored a local copy of a^^j), /S^j-^, and of 
the k-th element of vectors Vf^,, Vy_i). Then, the normalization step (14)-(15) is similar to the 
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Algorithm 1: Decentralized power method 
Input : Received signal vectors y[A;] G for k E {1, . . . , K}; number of iterations M; 

starting values f [k] V/c; averaging time t. 

Output: Eigenvalue estimates Ai[A;] \fk. 

1 for all nodes k in parallel do 

2 for iteration j = 1 to M do 

z[A;]^ = AC*vMK-.,W*-y[A:]^); 

4 Compute locally i\j-)[k] = ^z[k]^y[k]; 

5 end 

6 z[kf = kC%[k]{v,,,,[k]*-y[kf); 

7 d[k]=AC\[k]{\v,,,,[k]n 

8 Compute locally Ai[A;] = f ■ \\z[k]\\y d[k]; 

9 end 



denominator of (8), therefore it can be implemented by scalar AC and element-wise division. 
Based on the above observations, we can state the following result. 

Proposition 3: Given an ideal AC routine AC.(-), the LA (12)-(15) can be rewritten as 

= ^RN[AC;v (Diag(v*,) ■ Y) ] (33) 

wu) = ^diag{Y ■ [ACAr(Diag(v*,) ■ Y)]""} - a^,,Y,^, - A,)Vo_i) (34) 

/35^,,lx = i^-ACi(w*,0W(,,) (35) 

V(,+i) = ^u)//3u+i)- (36) 

Proof: The proof is a combination of the same steps already used in Prop. 1 and Prop. 2. 
More precisely, (33) follows from (12) using (26)-(30); (34) from (13) using (17); (35) from 
(14) using (31)-(32); and (36) is simply (15). □ 
The above result can be mapped to the decentralized algorithm reported in Alg. 2, where a 
realistic AC scheme AC* is adopted. 
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Algorithm 2: Decentralized Lanczos algorithm 



Input : Received signal vectors y[A;] G for A; G {1, ... , K}; number of iterations 
M <K; starting values v^,^[k] VA;, such that Ylk=i K'cdWP = 1' AdM = 
averaging time t. 

Output: Eigenvalue estimates {Ai[A;], . . . , Aj\f [A;]} V/c. 



1 for all nodes k in parallel do 



2 
3 
4 
5 
6 
7 
8 
9 
10 



for iteration j = 1 to M do 



Compute locally a^^) \k] 



El 

N 



Compute locally w^^^yk] = §z[k]"y[k] - a^^^[k] ■ v^^^[k] - (3y^[k] ■ v^^_^^[k]; 
b[k] = AC\[k]{\w,,M'y^ 

Compute locally /3(,+i)[A;] = ^/K ■ b[k] and v^,+^^[k] = w^,^[k]/ (3^,+^[k]; 



end 



Construct locally T[A;] from (11) using . . . ,a^M)[k], /3(2)[^], • • • , l3(M)[k]', 

Compute locally {\i[k], . . . /XM[k]} = eigenvalues of T[A;]; 



11 end 



V. Error and Complexity Analysis 

In this section we analyze the impact of non-ideal AC algorithms on the error and numerical 
complexity (especially in terms of network signaling) for the DPM and the DLA. 

A. DPM Error 

The DPM algorithm (Alg. 1) involves three sources of error due to AC. We define the first 
error term as 

E;^^.^^ 4 AC^v(Diag(v*_ J ■ Y) - lll^Diag(v* _ J ■ Y G C^^^, (37) 

which occurs at every iteration of (17), corresponding to line 3 of Alg. 1. Following the previously 
used convention, we denote by e|^^^[A;]-^ the A;-th row of E^Jf^- The second and third error terms 
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arise from (24), i.e., lines 6 and 7 of the algorithm, and are defined as 

4 AC^(Diag(v*,J ■ Y) - -^ll^Diag(v*,,,) • Y G C^x^, (38) 

e^^^ 4 ACUv*,,, © V(„,) - -^ll^(v;,,, © V(,,)) G C^. (39) 

Again, we refer to the A;-th row of E'^^^ as e^Jf^[k]'^ and to the A;-th element of e'^'^^ as e'^^^[A;]. 
Note that the three above defined errors are all instances of the general formula (4) applied with 
different inputs. To simplify the notation, we have dropped subscript t in the symbols of error 
variables. 

With regard to the DPM convergence, the most important term is clearly E^'^^, because this 
error is added at each iteration of the algorithm. The impact of E'^'^-^ on the evolution of PM 
vectors V(j, is expressed by the following result. 

Proposition 4: Given a non-ideal AC scheme that introduces an error term E™^ as defined 
in (37), the resulting DPM vector after M iterations can be written as 

M 

V,,, = R- v„ + 1 R-"Miag [Y (E-^) . (40) 

Proof: By applying the same steps as in the proof of Prop. 1 with the AC function defined 
in (4) with E^ replaced by E[^^^, we can write for any iteration j 

Va) = Rv,_, + ^diag[Y(Er^f)^]. (41) 

The above formula, applied recursively for M iterations, yields (40). □ 
The term R*^V(o) in (40) represents the ideal PM output, while the summation on the r.h.s. 
represents the error. For brevity, we define 

d,,4i-diag[Y(E-^)"], (42) 

which represents the error introduced by AC in a single iteration j. Its A;-th element (i.e., the 
component for node k) is d^j-j[k] = jj:{e^j^[k])^y[k]. 

Now, the convergence of the DPM depends on the relative magnitude of the error term 
JZjLi R-^^ ''d(j) compared to the desired term R^^V(„) as M — oo (recall that both terms are 
unnormalized). Let us assume that the spectral radius of R (i.e., Ai, since the eigenvalues are 
real and positive) is larger than 1."^ Then, the DPM error converges asymptotically to zero if 

''This assumption simplifies tlie mathematical analysis, and it is not a limitation in practice, because it is always possible to 
rescale the data samples Y such that Ai > 1. 
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the magnitude of the AC error vector grows slower than a certain rate. The result is expressed 
formally by the following proposition. 

Proposition 5: Let 6'(j, G [0,7r/2] be the angle between the true eigenvector and its DPM 
estimate V(j, after j iterations, defined by 

" -, 0<j<M, (43) 



cos Oi^j) = 



H ^U) 



'(J) I 



and assume cos^(o) 7^ 0, Ai > 1. 

Then, asymptotically in M, provided that ||d(i,^)||oo 



Ai 



A/+1 I max{A2,l} 



M 



we have 



lim I sin 6', 



(A/) I 



0. 



(44) 



Proof: Similarly as in [22, p. 406], we start by expressing vectors V(o) and d(j, in the 
eigenbasis (ui, . . . , uk), so that 



'(0) 



Oo,iUi + . . . + ao kUk, 



(45) 
(46) 



d(,) = Oj- lUi + . . . + ajj^UK, I < j < M. 

By assumption we have |ao,i| = cos 6^(0) 7^ 0. The DPM vector after M iterations can be then 
written as 



M 



M 

= ^aj,iAi^~-'ui + . . . ^ 



A4 



E\M~j 



j=0 



j=0 



and consequently 



sin^( 



M)| 



|ufV(jv,)| 



'(M)| 



1 - 



EM xi\ 



i=l l^j=o'^j,i^ 



M 



1=2 l^i=0^j,i'^ 



M 



M-j 
2 



i=l l^j=o'^j,i' 



M 



EK sr^M \n 
i=2 2^j=0 

2^j=o «j,i^i 



M 



< 



,M-j 



(47) 
(48) 

(49) 
(50) 

(51) 
(52) 
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By letting ai{M) = maxo<j<M {cLj^il and dividing numerator and denominator by , we have 

(52) < — 



EM ii 



E 



K 2M 
i=2 "i 



Er=o(A./Ai)^^Ar^ 



(53) 



(54) 



Y!f=o %,i'^i ' with Ai = 1/Ai e 



Let us first consider the denominator, which can be written as 
(0,1). We have to consider two cases, (i) Assume the series is absolutely convergent, so that 
livsiM^oo Ejlo I'^i.il ^ Moreover, since ao,i 7^ 0, we have limjv/^oo Ejlo I'^i.il ^ some 
a > 0. Now note that Ejio'^i.i-^i Z-transform of the sequence {aj_i}°^g at Ai. So, as 

Ai < 1, the Z-transform exists, and therefore the series converges to a value (3. (ii) Assume now 
that limM^oo Ejio I'^i.il ~ ^^^^ case, the radius of convergence of the power series in 

Ai is 0, which means that the series diverges for any value of Ai. Combining the two cases. 



we conclude that lim 



M-i^oo 



Ej=o ^ C)o], which means that in the worst case the 



denominator converges to a finite positive value. 
Consider now the numerator, which we can write as 



K 



A 



M 



^a,(M)^ 



i=2 



M 



j=0 



(55) 



Again, we have two cases, (i) If Aj > 1, then A^ > 1 for any < j < M, hence 

2 



M 
j=0 



M 



|(M + 1)(A,/Ai 



a/ 1 2 



and, recalling that Aj < Ai for any i>2 and applying L'Hopital's rule, we obtain limM-i>oo 
0. (ii) If Ai < 1, we have 



(56) 





M 


2 


M 


— 




< 


w 




i=o 




j=0 



|(M + i)Ar 



(57) 



whose limit is again 0. Combining these results yields 



K 



Am < 



i=2 



(M + i; 



/ max{A2, 1} 



AI 



(58) 
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Now suppose that max2<i<ft: aj(M) < 7(M). Then, the numerator converges to if 7(M) = 

I7TT ( max{A2,i} ) ^- '^^^ ^^'^^ ^^^^ max2<i<K ai{M) is equal (up to a constant) to ||d(„)||oo 
completes the proof. □ 
We next consider the impact of error terms E'^'^^ and e^'^^, which only concem the eigenvalue 
computation phase at the final iteration. Let enumi^] and eden[^] be the errors introduced by non- 
ideal AC in step 8 of Alg. 1, defined such that 

= V + en \k] • ^^^^ 

The following proposition provides the values of Cden [k] and Cnum [k] as a function of E^'^^ and 

gPM3_ 

Proposition 6: Given a non-ideal AC scheme that introduces error terms E'^'^^ and e'^'^^ 
defined in (38) and (39), the resulting errors in Ai[A;] are given by 

|2 



enu.M = ^[(e™^[A:])^Y%„, + v^„Y(e^^^[A:])*] +^||e^^^[fc]||^ (60) 
eden[A;] = i^-e™^[fc]. (61) 
Proof: Using the results of Prop. 2 and introducing non-ideal AC, we can write 



^iW = l^- "\ I. , .PM3 i" (62) 



ilK)Y|r + f [(e™2[A:])^Y%,,, + v^,,Y(e™2[A:])*] + f ||ePM2[fc] ||2 



(63) 



Since -i:||v^^jY||^ = v^^jRv(jv^), we note that (63) is equivalent to (59), with error expressions 
enum[^] and eden[^] given, respectively, by (60) and (61). □ 
By simple algebraic manipulations, and letting = (v^^jRV(jv^))/||V(a/)P be the estimate 
of Ai at the A'/-th DPM iteration without eigenvalue computation errors, we can write 

^^[^] = TT p ■ ^1 + II 'Zf^ i,r (64) 

1 + edenFj/||V(M)r \\^{M)\\ + CdenFj 

The above expression shows immediately that the eigenvalue computation error becomes asymp- 
totically negligible, as ||v(a/)|P — )■ oo as M — )■ oo. 
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B. DLA Error 

The DLA involves two error terms due to AC, in lines 3 and 6 of Alg. 2. The first error arises 
from in (33)-(34) and is defined as 

4 AC5v(Diag(v*,) ■ Y) - -ill^Diag(v*,) ■ Y G C^^^. (65) 

The second error originates from (35) and is defined as 

4 AC^w*, w,,) - ^11^ (w*,, w,,) G C^. (66) 



As usual, the /c-th row of is denoted by e)^^'^[kY , and the k-\h element of e[-j^^ by e[-^!^^[/c]. 
Both error terms occur at every iteration j of the algorithm. We are now interested in the impact 
of such errors on W(j). First of all, we notice that the /c-th component of vector W(j) in the ideal 
LA (13) can be written as 

<M ^ |y[/c]^Y^Vo) - (vJ,Rv,,) ■ v,,M - ||W(,_,,|| ■ v,,_.,[kl (67) 

by exploiting the expressions of (12) and (14). We then define the error on to be 



.J [k] = [k] - wjj) [k] . (68) 

A closed-form expression of e^^(^j)[k] as a function of E|-^'^^ and e^-^'^^ on e^(j)[/i;] is given by 
the following Proposition. Later, we provide an interpretation of this error expression, and we 
discuss its impact on the estimation of eigenvalues. 

Proposition 7: Given a non-ideal AC scheme that introduces error terms E^-^'^"^ and e\^^^^ defined 
respectively in (65) and (66), the resulting error in if(j)[A;] is given by 

e.U)[k] =^ie\,';'[kry[k] - I [(eLAi[A;])^Y^v,,,, + vJ,,Y(e;;^^[A:])*] ■ v,,,[k] 

- ^H!^"[kW ■ W - ^^^^0?'^ ■ ^0-) W + oie\^'[k]). (69) 

Proof: The LA iteration (13), as well as the decentralized version (34), consists of three 
additive terms, therefore the error (68) can be expressed as the sum of three separate terms: 

e.Jk] = e'^Jk] + e^l[k]+e'^^;i[k]. (70) 

The first term originates from the computation of ^diagjY ■ [ACjv(Diag(v* j) ■ Y)]^} in (34), 
which is identical to the PM iteration. Thus, the first error term can be expressed using (41) 
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with E^Jf^ replaced by E[-j^^. The global error for all nodes is ;^diag[Y(E[-j^^)^], and its k-th 
element is 

e^lfi] = j^{e\,^'[kry[k]. (71) 

The second term arises from computation of which is done using the same vector z[k] 

obtained via AC in line 3. Therefore the error on depends again on E^-^"^-^. Since 

is calculated as the square norm of z[k] (line 4), the structure of the error is the same as that 
of enum[^] for the DPM (60). By replacing E'^'^^ by Ej-^'^^ in (60) and multiplying by i\j)[k], we 
obtain the second part as 

eZik] = f [(e;;^^[A;])^Y%„, + v;i,^Y(e;;f [fc])*] ■ v,,[k] + ^He^^f ■ (72) 

The third part contains the error due to computation of (3(^j)[k], i.e., the norm of W(j_i). The 
error arises from the AC algorithm used in line 6 when computing ||W(j_i)|p; a nonlinearity is 
then introduced by the square root. Using a first-order Taylor expansion (under the assumption 
that e[-^'^^[A;] ^ ||), we can write 

/3o)W = \/l|wo_.,P + i^eLfW (73) 

(74) 

e\^'[k] + o{e\^'[k])^ (75) 
'■[k]+o{e\;^'[k]), (76) 

^11^0-1)11 

hence the third error term is given by 

The resulting error (69) is obtained by summation of (71), (72), and (77). □ 
The error expressed by Prop. 7 is then propagated to the next iteration j + 1 through another 
nonlinear step (line 7). Using the same procedure of Eqs. (73)-(76), the relation between fy+i)[A;] 
and can be written as 

In summary, (69) expresses the error on W(^) given V(j, and v^-i), and (78) expresses the error 
on V(j+i) given w^^y In principle, a closed-form error update rule could be derived from these 
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expressions (as it was done for the DPM), but the resuking formula would be too complicated 
to provide additional insight. We conclude the DLA analysis with two remarks. 

1) Following the same line of reasoning of [26], the error term in (69) represents a loss of 
orthogonality of the vectors v^^j, while (78) results in a loss of unit norm. As documented in the 
literature, the loss of orthogonality leads to the appearance of so-called spurious eigenvalues at 
certain iterations (typically spurious eigenvalues are duplicates of existing eigenvalues already 
computed in the previous iterations). Some heuristics for the identification of spurious eigenvalues 
exist, such as the CuUum-Willoughby method [40], which compares the eigenvalues of T (11) 
with those of the same matrix without the first row and column. Alternatively, spurious values 
can be detected by comparing the eigenvalues of T at iterations j and j — 1 . In our application, 
another simple criterion can be derived by observing that the covariance matrix R has exactly 
min{i^, A^} non-zero eigenvalues, hence for iterations j > min{i^, A^} any new non-zero value 
is automatically identified as spurious. We remark that all criteria for detection of spurious 
eigenvalues are compatible with the DLA, as they involve local post-processing of the output at 
individual nodes (matrices T[k] as defined in line 9 of Alg. 2). 

2) The terms ^ A !^^[k] ■ v,^_u[k] in (69) and ^^ej-^^ (yg) may be critical for 
the convergence of the algorithm in the event that ||W(j)|| ^ (or, equivalently, ^ 0) at 
some iteration j. In the ideal LA, quoting [22], "a zero (3^j) is a welcome event in that is signals 
the computation of an exact invariant subspace. However, an exact zero or even a small 

is a rarity in practice." This event is even more unlikely to happen in the case of the DLA: 
simulation results confirm that the values of (3^^^ are always far from zero. As a result, cases 
of divergence of the DLA were never observed. Nevertheless, the DLA turns out to be more 
sensitive to numerical problems than the DPM, as shown in Section VII. 

C. Complexity 

For both the DPM and the DLA, complexity mainly arises from the repeated use of AC 
routines, resulting in communication overhead, time delays, and possible synchronization issues. 
We next compare the complexity of the two algorithms in terms of the following parameters: (i) 
number of calls to functions ACn and ACi; (ii) total number of "information units" exchanged over 
the wireless channel by one node with degree (number of neighbors) d, where an information 
unit is defined as the number of bits used to encode a scalar; (Hi) number of algorithmic steps. 
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Number of 


Number of 


Number of information 


Number of required 




ACat 


ACi 


units exchanged per node 


time periods 


DPM 


M + 1 


1 


I{MN + N + l)d 


AI + 2 


DLA 


M 


M 


I{MN + AI)d 


2M 



TABLE II 

Numerical complexity of DPM and DLA. Legend: TV = number of samples, M = number of DPM/DLA 

ITERATIONS, I = NUMBER OF AC ITERATIONS, d = NODE DEGREE. 



defined as individual tasks that must be performed in a sequential way due to input-output 
dependency (i.e., step n cannot start until step n — 1 is completed). 

The above performance figures are reported in Table II for the two algorithms. For simplicity 
it is assumed that all instances of AC have the same number of iterations, /. It can be observed 
that the total number of calls to consensus routines is slightly higher for the DLA (assuming 
M > 2): the DPM requires one vector- AC per iteration (line 3) in addition to another vector- AC 
and a scalar-AC for eigenvalue computation (lines 6-7), whereas the DLA uses one vector-AC 
and one scalar-AC in each iteration (lines 3-6). The amount of information sent over the air is 
nearly the same, as shown in the third column of the table. However, the DLA is significantly 
more complex than the DPM in terms of the number of algorithmic steps: in the case of the 
DLA, each iteration consists of two sequential steps (vector computation and normalization) that 
cannot be parallelized, hence the total number of steps for the DLA is nearly twice that of the 
DPM. 

The above results suggest that the values of M, N, and / should be kept as small as possible 
in order to reduce complexity and improve the reactivity of the algorithms, especially when used 
in real-time detection applications. In particular, the number of samples can be very small if 
the number of nodes K is large enough (large networks are the natural scenario of application 
for decentralized algorithms). For example, in a network of 40 nodes, 10 samples per node are 
sufficient to detect a signal with SNR of 7dB with high probability (see Fig. 5(b) in Section 
VII). 

We now briefly analyze the computational complexity for individual nodes. For the DPM, 
the dominant factor is the vector product of line 4. Since the vector size is and the product 
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Name 


Test statistic 


Application scenario 


Ref. 


Roy's test (RT) 


Tr = Ai/a2 


P = 1, known 


[30], [32], [37] 


GLR test (GT) 




P = 1, unknown 


[34]-[36] 


Sphericity test (ST) 




P >1, finite SNR 


[34], [35], [38] 


John's test (JT) 




P > 1, low SNR 


[38] 



TABLE III 

Eigenvalue-based tests for multi-sensor signal detection. 



is computed at every iteration, the computational complexity per node scales as 0{MN). In 
the DLA, every iteration involves computing the norm of a vector of size N (line 4) and a 
vector product of the same size (line 5), hence the computational complexity per node scales 
as 0{MN) as well. In addition, the DLA involves local calculation of the eigenvalues of the 
tridiagonal matrix T[A;] (line 10), which has complexity 0{K'P). The dominant term between 
0{MN) and 0{KP) depends on the relative values of M and N. 

VI. Application to Spectrum Sensing 

We consider a distributed spectrum sensing scenario, where multiple sensor nodes cooperate 
to detect the presence of a primary signal in a given frequency band. A decision about signal 
presence or absence is made upon receiving N signal samples at each sensor. The main challenge 
for cognitive networks is to achieve reliable signal detection with a limited number of samples. 
IVIathematically, the problem is a binary hypothesis test between a "signal-plus-noise hypothesis" 
(Hi) and a "noise-only" hypothesis (Hq) based on the received samples Y. The network is 
assumed to operate under homogeneous conditions, i.e., the same hypothesis {Hq or Hi) holds 
for all sensors during the entire sensing period. Given the model already introduced in Section 
II and assuming zero-mean complex Gaussian noise, the signal vector received at a given time 
instant n at the K sensors can be written under Ho as 

y(n)|^„ = r7(n), (79) 

where 77(77.) ~ A/'c(0, cr^I). Under Hi, the received vector is 

yHIh, =Hs(r7) + 77(77), (80) 
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where s{n) ~ A/c(0, Diag{[aj, . . . , cxp])) G is a vector of signal samples transmitted by P 
sources (primary users), modeled as zero-mean Gaussian mutually uncorrelated random variables, 
and H = [hi,..., hp] E (^kxp ^ complex channel matrix whose columns represent the 
channels between the P signal sources and the K sensors. The channel coefficients are assumed 
to be unknown but constant during the sensing period^. The SNR for the i-th signal source is 
defined (under "Hi) as pi = ||hj|p(Tj^/(T^. Let T(Y) be a generic test statistic for signal detection, 
computed from the signal samples, and let be the associated decision threshold, such that the 
detector decides for T-Li if T(Y) > ^9 and for "Hq otherwise. Then, false-alarm and detection 
probabilities are defined, respectively, as Pfa = Pr[T(Y) > i^lKo] and Pj = Pr[T(Y) > i^lKi]. 
Detectors are typically calibrated so as to achieve a fixed false-alarm rate Pfa = a, i.e., the 
threshold is set as 

Signal detection in the above-described setting has been extensively studied in the cognitive 
radio literature, e.g., [30]-[38], leading to the derivation of several possible test statistics T(Y). 
Most of such statistics are functions of the eigenvalues of the sample covariance matrix (some 
examples are reported in Table III) and, as such, they can be computed in a decentralized network 
by applying the DPM or the DLA proposed in this paper. The problem reduces to computing, 
for each sensor k, a local version of the adopted test statistic, say Ti[k] with i E {R, G,S,J}. 
This is obtained directly as a function of the local eigenvalue estimates \j[k] (with j = I for 
the DPM, which is sufficient to compute the RT statistic; and I < j < K for the DLA). Then, 
each node tests the local statistic Ti[k] against the predefined threshold '&{a) and decides for "Hq 
or "Hi according to the rule 

f,[k]'^^{a). (81) 

-Ho 

In order to average out the numerical errors introduced by non-ideal AC at different nodes^, a 
final round of AC can be executed taking as inputs the local statistics Ti[k]. The result 

f:[k]=AC{[k]{f,[k]) (82) 

shall be used in (81) instead of Ti[k]. 

^Assuming the sensing period to be shorter than the coherence time of the channel is reasonable in multi-sensor spectrum 
sensing applications, where accurate detection is achieved already at low sample size. 

^Numerical errors may result in the problem of different decisions at different nodes if Ti[k] is very close to ^{a). 
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Fig. 1. Network topology. 



The popular cooperative energy detector, i.e., the test based on the (possibly weighted) sum of 
the received signal energies at different sensors, also admits a natural decentralized implemen- 
tation via AC algorithms. This problem was investigated in [39]. Decentralized energy detection 
is computationally simpler than eigenvalue-based techniques, but clearly inherits the well-known 
shortcomings of energy detection (suboptimality in multi- sensor settings and sensitivity to noise 
uncertainty). 

In some applications, the goal is not only to discriminate between Hq and Hi, but to estimate 
the number of signals P. This problem can be addressed again by eigenvalue-based estimators, 
such as the well-known Wax and Kailath's information-theoretic criteria [41], or the recent 
random matrix theory (RMT) estimator proposed in [42]. In all cases, the estimated number of 
signals is estimated as 

P = arg min T(g;Ai,...,Ai^), (83) 

0<g<il 

where T(-) is a function of the eigenvalues and, once again, can be computed locally from the 
DLA estimates Ai, . . . , Ai^. 

VII. Simulation Results 

For the purpose of numerical evaluation of the proposed algorithms, we consider a randomly 
generated network consisting of K = AO nodes, as depicted in Fig. 1. Edges in the graph represent 
communication links between pairs of nodes. The received signal samples Y are randomly 
generated at every Monte Carlo iteration according to the signal-plus-noise model described in 
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5 10 15 20 25 30 

AC iterations 



Fig. 2. Comparison of different AC algorithms. 

Section VI (case "Hi), with the following parameters: = 10 samples per node, noise variance 
0"^ = 1, one Gaussian signal source (P = 1) with SNR p = 5dB. 

We first turn our attention to the convergence of the largest eigenvalue, expressed in terms 
of the mean square error (MSE) of Ai with respect to Ai. In practice, the MSE is calculated 
as an average over 3000 Monte Carlo simulations. In Fig. 2 we compare the performance of 
different AC algorithms when applied in the DPM (Alg. 1). The number of DPM iterations is 
fixed to M = 20, while the number of AC iterations (/) varies from 5 to 30 (x-axis in the 
plot). We consider four alternative AC schemes, namely two versions of traditional AC [8] - 
with Metropolis weights (a heuristic based on the graph Laplacian) and with optimal symmetric 
weights (resulting from the solution of a convex optimization problem) - and two versions 
of AC with Chebyshev acceleration [19] - assuming first a fixed network topology and then 
one with 3% link failure probability. The lowest achievable bound is given by the ideal PM 
(i.e., a DPM with exact AC). The Chebyshev acceleration algorithm turns out to significantly 
outperform traditional AC even with optimized weights (under deterministic settings). For this 
reason we adopt Chebyshev-accelerated AC as the standard consensus scheme in all the following 
simulations. 

In Fig. 3 we compare the DPM against the DLA. We observe that the DLA exhibits a faster 
convergence rate, but is more sensitive than the DPM to numerical errors introduced by imperfect 
consensus: the DLA error converges to for J = 15 AC iterations, but not for / = 10, whereas 
the performance of the DPM is practically identical to that of the ideal PM already for / = 10 



March 29, 2013 



DRAFT 



25 




2 4 6 8 10 12 14 16 

Algorithm iterations 



Fig. 3. Convergence of Ai. 

(and even for smaller values). In other words, the performance of the DLA compared to the 
DPM is better in absolute terms (especially with few algorithm iterations), but worse in relative 
terms. This behavior is consistent with the analysis presented in Section V and can be understood 
intuitively by noticing that each step of the DLA uses AC twice, in contrast to the DPM where 
AC is used only once per iteration. 

We next evaluate the performance of the DLA when estimating multiple eigenvalues. The 
results, shown in Fig. 4, indicate that the higher is the order of eigenvalues to be estimated, 
the more iterations are needed. The reason is that, by definition of the LA, the j-th eigenvalue 
cannot be estimated until the j-th iteration; in addition, numerical errors sometimes cause the 
appearance of spurious eigenvalues (see Section V-B) which, even when correctly identified, 
introduce a one-step delay in the algorithm. One way to mitigate the numerical problems that 
affect the estimation of high-order eigenvalues is to improve the precision of the AC routine 
by tuning the parameter /. For example, according to our simulations, 20 AC iterations are 
sufficient to achieve an accurate estimation of the Ai and A3 (Fig. 4(a)), whereas 30 iterations 
are necessary for A5 and Ag (Fig. 4(b)). 

We now illustrate an application of the proposed DPM and DLA for spectrum sensing in a 
distributed cognitive radio network. We assume again the same network topology of Fig. 1, with 
K = AO and = 10. According to the model introduced in Section VL the samples under "Ho 
are Gaussian distributed with variance cr^ = 1, while under "Hi we have one Gaussian signal 
component (P = 1) with SNR p = 7dB. We test the performance of two signal detectors: the RT 
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Lanczos iterations (M) Lanczos iterations (1^) 

(a) / = 20 (b) 1 = 30 

Fig. 4. Convergence of multiple eigenvalues. 

and the GT, defined respectively as Tr and Tq in Tab. III. The RT is a test of significance of the 
largest eigenvalue (Ai) alone, hence it can be implemented in a decentralized setting using either 
the DPM or the DLA. The GT, in contrast, involves all eigenvalues and requires the use of the 
DLA. In Fig. 5 we compare: (i) the ideal performance of RT and GT using the eigenvalues of 
R computed exactly by a fusion center with perfect communication links; (ii) the performance 
achieved when the eigenvalues are computed by the PM or LA (assuming ideal AC), respectively 
for M = 5 and M = 10; and (Hi) the performance achieved when the eigenvalues are computed 
by the DPM or DLA, with / = 30 iterations. The results show that, after M = 5 algorithm 
iterations (Fig. 5(a)), the RT using LA or DLA already attains nearly-ideal performance, while 
the RT with DPM is slightly suboptimal; for the GT, on the other hand, 5 iterations are not 
enough to reach convergence of all eigenvalues (note that the performance gap is due to the LA 
itself and not to the decentralized version). After M = 10 iterations (Fig. 5(b)), all versions of 
the RT converge to the ideal RT bound, and the gap of LA- and DLA-GT compared to the ideal 
GT is dramatically reduced. 

VIII. Conclusion 

In this paper we have proposed and analyzed two general-purpose algorithms that can be used 
for computing sample co variance eigenvalues in distributed wireless networks. As an application, 
we have considered spectrum sensing in cognitive networks, and we have shown that numerous 
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Fig. 5. Signal detection: ROC curves. 
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eigenvalue-based tests for single-signal and multiple-signal detection can be implemented in 
a decentralized setting by using the proposed algorithms. Such decentralized signal detection 
techniques enable sensor nodes to compute global test statistics locally, thereby performing 
hypothesis tests without relying on a fusion center. Decentralized approaches also provide 
additional robustness to node failures or Byzantine attacks. 

The two proposed algorithms - the DPM and the DLA - have different strengths and draw- 
backs. The DPM is less complex, more robust to numerical problems, and provably convergent 
even in the presence of non-ideal AC (under the conditions of Proposition 5); however, it can 
estimate only the largest eigenvalue, and under ideal assumptions its convergence rate is slower 
than that of the DLA. On the other hand, the DLA is able to estimate all eigenvalues (albeit 
with increasing complexity with the number of eigenvalues) and, even in the presence of AC 
errors, provides faster initial convergence than the DPM; however, it is more complex and more 
sensitive to numerical errors, and requires some post-processing in order to remove "spurious" 
eigenvalues. 

Evaluated through simulations, both algorithms exhibit good performance in practical con- 
ditions (non-ideal AC algorithms, small number of samples) already after few iterations. In 
particular, convergence is very fast in the case of the largest eigenvalue, which results in 
high-performing distributed signal detectors based on the largest eigenvalue (referred to as RT- 
DPM and RT-DLA in Fig. 5). Other possible fields of applications of the proposed algorithms 
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are distributed anomaly detection in wireless sensor networks and signal feature estimation in 
distributed antenna arrays. 
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